rm(list = ls())
#setwd("YOUR-WORKING-DIRECTORY")

library("list")
library("rr")

## Set seed
set.seed(4444)

## Load data
load("dir-data.RData")
load("list-data.RData")
load("rr-data.RData")

## Common covariate formula
# Include these variables:
# age, age squared, gender, PID, education
cov.formula <- formula(~ ageyrs.new + ageyrs2.new + male + pidREP + pidIND + eduHIGHER)

#################
##   DIRECT    ##
#################

dir.logit <- glm(update.formula(cov.formula, dir ~.), family=binomial(link=logit), data=dir)
summary(dir.logit)

vars <- c("male", "pidREP", "eduHIGHER")
model.matrix <- data.frame(model.matrix(dir.logit))
predict <- NULL
tab <- matrix(NA, 3, 6)
rownames(tab) <- vars
colnames(tab) <- c("est.1", "low.1", "hi.1", "est.0", "low.0", "hi.0")
j <- 0

for (i in vars){
  j <- j+1
  new.model.matrix <- subset(model.matrix, get(i)==1)
  predict <- predict(dir.logit, newdata=new.model.matrix, type="response", se.fit=TRUE)
  tab[j,1] <- 1-mean(predict$fit) # predicted probability of "no" vote
  tab[j,2] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
  tab[j,3] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)
  new.model.matrix <- subset(model.matrix, get(i)==0)
  predict <- predict(dir.logit, newdata=new.model.matrix, type="response", se.fit=TRUE)
  tab[j,4] <- 1-mean(predict$fit) # predicted probability of "no" vote
  tab[j,5] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
  tab[j,6] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)
}

# Replace prediction for non-republicans with prediction for democrats
new.model.matrix <- subset(model.matrix, pidREP==0 & pidIND==0)
predict <- predict(dir.logit, newdata=new.model.matrix, type="response", se.fit=TRUE)
tab[2,4] <- 1-mean(predict$fit) # predicted probability of "no" vote
tab[2,5] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
tab[2,6] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)

dir.tab <- tab
rownames(dir.tab) <- c("male", "pidREP", "eduHIGHER")



###############
##   LIST    ##
###############

ml <- ictreg(list.y ~ ageyrs.new + ageyrs2.new + male + pidREP + pidIND + eduHIGHER, 
             data = list, 
             treat = "treat", 
             J = 4, 
             method = "ml",
             constrained = TRUE,
             fit.start = "lm")
summary(ml)

vars <- c("male", "pidREP", "eduHIGHER")
model.matrix <- data.frame(ml$x) 
predict <- NULL
tab <- matrix(NA, 3, 6)
rownames(tab) <- vars
colnames(tab) <- c("est.1", "low.1", "hi.1", "est.0", "low.0", "hi.0")
j <- 0

for (i in vars){
  j <- j+1
  new.model.matrix <- subset(model.matrix, get(i)==1)
  predict <- predict(ml, newdata=new.model.matrix, avg=T, se.fit=T)
  tab[j,1] <- 1-mean(predict$fit) # predicted probability of "no" vote
  tab[j,2] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
  tab[j,3] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)
  new.model.matrix <- subset(model.matrix, get(i)==0)
  predict <- predict(ml, newdata=new.model.matrix, avg=T, se.fit=T)
  tab[j,4] <- 1-mean(predict$fit) # predicted probability of "no" vote
  tab[j,5] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
  tab[j,6] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)
}

# Replace prediction for non-republicans with prediction for democrats
new.model.matrix <- subset(model.matrix, pidREP==0 & pidIND==0)
predict <- predict(ml, newdata=new.model.matrix, avg=T, se.fit=T)
tab[2,4] <- 1-mean(predict$fit) # predicted probability of "no" vote
tab[2,5] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
tab[2,6] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)

list.tab <- tab
rownames(list.tab) <- c("male", "pidREP", "eduHIGHER")


#### 
load("ms-g11-county-vote-data.Rdata")

# Select sampled counties from voting data
vote <- vote[vote$insample == 1, ]

# Create numeric group indicators
list$insample <- rep(1,nrow(list))
list$numeric.county <- as.numeric(list$county)

# define auxiliary moments (in terms of 'yes vote)
# vote return for all 19 counties together
aux.state <- with(vote, sum(YES_TOT)/sum(size))
# vote return for each of 19 counties (vector of county returns in alphabetical order)
aux.counties <- vote$yes26

## ESTIMATES WITH COUNTY-LEVEL CONSTRAINTS
# County constraints
# setting aux moment
h <- aux.counties
names(h) <- 1:19
# setting group
group <- list$numeric.county


aux <- ictreg(list.y ~ ageyrs.new + ageyrs2.new + male + pidREP + pidIND + eduHIGHER, 
             data = list, 
             treat = "treat", 
             J = 4, 
             method = "nls", 
             h = h, 
             group = group,
             fit.start = "lm")
summary(aux)

vars <- c("male", "pidREP", "eduHIGHER")
model.matrix <- data.frame(ml$x) 
predict <- NULL
tab <- matrix(NA, 3, 6)
rownames(tab) <- vars
colnames(tab) <- c("est.1", "low.1", "hi.1", "est.0", "low.0", "hi.0")
j <- 0

for (i in vars){
  j <- j+1
  new.model.matrix <- subset(model.matrix, get(i)==1)
  predict <- predict(aux, newdata=new.model.matrix, avg=T, se.fit=T)
  tab[j,1] <- 1-mean(predict$fit) # predicted probability of "no" vote
  tab[j,2] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
  tab[j,3] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)
  new.model.matrix <- subset(model.matrix, get(i)==0)
  predict <- predict(aux, newdata=new.model.matrix, avg=T, se.fit=T)
  tab[j,4] <- 1-mean(predict$fit) # predicted probability of "no" vote
  tab[j,5] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
  tab[j,6] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)
}

# Replace prediction for non-republicans with prediction for democrats
new.model.matrix <- subset(model.matrix, pidREP==0 & pidIND==0)
predict <- predict(aux, newdata=new.model.matrix, avg=T, se.fit=T)
tab[2,4] <- 1-mean(predict$fit) # predicted probability of "no" vote
tab[2,5] <- (1-mean(predict$fit))-qnorm(.975)*mean(predict$se.fit)
tab[2,6] <- (1-mean(predict$fit))+qnorm(.975)*mean(predict$se.fit)

aux.list.tab <- tab
rownames(aux.list.tab) <- c("male", "pidREP", "eduHIGHER")



##############################
##   RANDOMIZED RESPONSE    ##
##############################

fit <- glm(update.formula(cov.formula, rr.dum ~.), data=rr, family=binomial)
start <- fit$coefficients
rr.mod <- rrreg(rr.dum ~ageyrs.new + ageyrs2.new + male + pidREP + pidIND + eduHIGHER, data = rr, p = 1/2, p0 = 0, p1 = 1/2, verbose=TRUE, start = start, maxIter = 10000, design = "forced-known")             
summary(rr.mod)

vars <- c("male", "pidREP", "eduHIGHER")
model.matrix <- data.frame(rr.mod$x) 
predict <- NULL
tab <- matrix(NA, 3, 6)
rownames(tab) <- vars
colnames(tab) <- c("est.1", "low.1", "hi.1", "est.0", "low.0", "hi.0")
j <- 0

for (i in vars){
  j <- j+1
  new.model.matrix <- subset(model.matrix, get(i)==1)
  predict <- predict(rr.mod, newdata=new.model.matrix, given.y=F, avg=T, quasi.bayes=T)
  tab[j,1] <- 1-predict$est # predicted probability of "no" vote
  tab[j,2] <- (1-predict$est)-qnorm(.975)*predict$se
  tab[j,3] <- (1-predict$est)+qnorm(.975)*predict$se
  new.model.matrix <- subset(model.matrix, get(i)==0)
  predict <- predict(rr.mod, newdata=new.model.matrix, given.y=F, avg=T, quasi.bayes=T)
  tab[j,4] <- 1-predict$est # predicted probability of "no" vote
  tab[j,5] <- (1-predict$est)-qnorm(.975)*predict$se
  tab[j,6] <- (1-predict$est)+qnorm(.975)*predict$se
}

# Replace prediction for non-republicans with prediction for democrats
new.model.matrix <- subset(model.matrix, pidREP==0 & pidIND==0)
predict <- predict(rr.mod, newdata=new.model.matrix, given.y=F, avg=T, quasi.bayes=T)
tab[2,4] <- 1-predict$est # predicted probability of "no" vote
tab[2,5] <- (1-predict$est)-qnorm(.975)*predict$se
tab[2,6] <- (1-predict$est)+qnorm(.975)*predict$se

rr.tab <- tab
rownames(rr.tab) <- c("male", "pidREP", "eduHIGHER")


# Create numeric group indicators
rr$insample <- rep(1,nrow(rr))
rr$numeric.county <- as.numeric(rr$county)
group <- rr$numeric.county[as.numeric(rownames(model.matrix))]

rr.aux <- rrreg(rr.dum ~ageyrs.new + ageyrs2.new + male + pidREP + pidIND + eduHIGHER, data = rr, p = 1/2, p0 = 0, p1 = 1/2, verbose=TRUE, start = coef(rr.mod), maxIter = 10000, design = "forced-known", 
  h = h, group = as.character(group))             
summary(rr.aux)

j <- 0
for (i in vars){
  j <- j+1
  new.model.matrix <- subset(model.matrix, get(i)==1)
  predict <- predict(rr.aux, newdata=new.model.matrix, given.y=F, avg=T, quasi.bayes=T)
  tab[j,1] <- 1-predict$est # predicted probability of "no" vote
  tab[j,2] <- (1-predict$est)-qnorm(.975)*predict$se
  tab[j,3] <- (1-predict$est)+qnorm(.975)*predict$se
  new.model.matrix <- subset(model.matrix, get(i)==0)
  predict <- predict(rr.aux, newdata=new.model.matrix, given.y=F, avg=T, quasi.bayes=T)
  tab[j,4] <- 1-predict$est # predicted probability of "no" vote
  tab[j,5] <- (1-predict$est)-qnorm(.975)*predict$se
  tab[j,6] <- (1-predict$est)+qnorm(.975)*predict$se
}

# Replace prediction for non-republicans with prediction for democrats
new.model.matrix <- subset(model.matrix, pidREP==0 & pidIND==0)
predict <- predict(rr.aux, newdata=new.model.matrix, given.y=F, avg=T, quasi.bayes=T)
tab[2,4] <- 1-predict$est # predicted probability of "no" vote
tab[2,5] <- (1-predict$est)-qnorm(.975)*predict$se
tab[2,6] <- (1-predict$est)+qnorm(.975)*predict$se

rr.aux.tab <- tab
rownames(rr.aux.tab) <- c("male", "pidREP", "eduHIGHER")

indiv.tab <- rbind(list.tab, aux.list.tab, rr.tab, rr.aux.tab)

#Set point sz/line to be used
x <- .8
a <- -2.5

pdf("figure-2.pdf", height = 4, width = 8.5) #open pdf device
par(mfrow=c(1,1), las=1, mar=c(0,2,2,1), oma=c(1,3,1,0), cex=1) 

x.coords <- c(c(1.33, 1.66, 1.99, 2.33), # Grouped by variable level
              c(3.33, 3.66, 3.99, 4.33),
              c(6.33, 6.66, 6.99, 7.33),
              c(8.33, 8.66, 8.99, 9.33),
              c(11.33, 11.66, 11.99, 12.33),
              c(13.33, 13.66, 13.99, 14.33))

plot(0:23, seq(from = -2, to = 1.5, length = 24), type='n', xlim=c(1,14.66), ylim=c(0,1), axes=FALSE,
     main="", ylab = "")
axis(2, at=seq(0, 1, .1), las=1, cex.lab=1.2)
par(las=0)
mtext("Estimated proportion of 'no' votes on Personhood", outer=TRUE, side=2, cex=1, at=0.5, line=1.5, padj=.1)

est <- c(indiv.tab[seq(1,12,3), 1],
         indiv.tab[seq(1,12,3), 4],
         indiv.tab[seq(2,12,3), 1],
         indiv.tab[seq(2,12,3), 4],
         indiv.tab[seq(3,12,3), 4],
         indiv.tab[seq(3,12,3), 1])

low <- c(indiv.tab[seq(1,12,3), 2],
         indiv.tab[seq(1,12,3), 5],
         indiv.tab[seq(2,12,3), 2],
         indiv.tab[seq(2,12,3), 5],
         indiv.tab[seq(3,12,3), 5],
         indiv.tab[seq(3,12,3), 2])

up <- c(indiv.tab[seq(1,12,3), 3],
        indiv.tab[seq(1,12,3), 6],
        indiv.tab[seq(2,12,3), 3],
        indiv.tab[seq(2,12,3), 6],
        indiv.tab[seq(3,12,3), 6],
        indiv.tab[seq(3,12,3), 3])

segments(x.coords, low, x.coords, up)
points(x.coords, est, pch=(rep(c(1, 16, 0, 15), 6)), col="black", bg="black", cex=1)
abline(h=0.653, lty="dashed")

mtext("Gender", outer=TRUE, cex=1.3, at=0.20, line=-1, font=2)
mtext("Party", outer=TRUE, cex=1.3, at=0.5, line=-1, font=2)
mtext("Education", outer=TRUE, cex=1.3, at=0.835, line=-1, font=2)

mtext("Male", outer=TRUE, cex=1, at=0.13, line=a, font=2)
mtext("Female", outer=TRUE, cex=1, at=0.26, line=a, font=2)

mtext("Republican", outer=TRUE, cex=1, at=0.44, line=a, font=2)
mtext("Democrat", outer=TRUE, cex=1, at=0.57, line=a, font=2)

mtext("No Higher", outer=TRUE, cex=1, at=0.76, line=a, font=2)
mtext("Higher", outer=TRUE, cex=1, at=0.89, line=a, font=2)

legend("bottomright", pch = c(1, 16, 0, 15), 
  legend = c("List Experiment", "List Experiment with Auxiliary Information", 
    "Randomized Response", "Randomized Response with Auxiliary Information"), 
  cex = x, 
  bty = "n")
dev.off()


## #######
## APPENDIX FIGURE 8
## #######

indiv.tab <- rbind(dir.tab, list.tab, aux.list.tab, rr.tab, rr.aux.tab)

#Set point sz/line to be used
x <- .8
a <- -2.5

pdf("appendix-figure-8.pdf", height = 4, width = 8.5) #open pdf device
par(mfrow=c(1,1), las=1, mar=c(0,2,2,1), oma=c(1,3,1,0), cex=1) 

x.coords <- c(c(1, 1.33, 1.66, 1.99, 2.33), # Grouped by variable level
              c(3, 3.33, 3.66, 3.99, 4.33) + 0.33,
              c(6, 6.33, 6.66, 6.99, 7.33),
              c(8, 8.33, 8.66, 8.99, 9.33) + 0.33,
              c(11, 11.33, 11.66, 11.99, 12.33),
              c(13, 13.33, 13.66, 13.99, 14.33) + 0.33)

plot(0:29, seq(from = -2, to = 1.5, length = 30), type='n', xlim=c(1,14.66), ylim=c(0,1), axes=FALSE,
     main="", ylab = "")
axis(2, at=seq(0, 1, .1), las=1, cex.lab=1.2)
par(las=0)
mtext("Estimated proportion of 'no' votes on Personhood", outer=TRUE, side=2, cex=1, at=0.5, line=1.5, padj=.1)

est <- c(indiv.tab[seq(1,15,3), 1],
         indiv.tab[seq(1,15,3), 4],
         indiv.tab[seq(2,15,3), 1],
         indiv.tab[seq(2,15,3), 4],
         indiv.tab[seq(3,15,3), 4],
         indiv.tab[seq(3,15,3), 1])

low <- c(indiv.tab[seq(1,15,3), 2],
         indiv.tab[seq(1,15,3), 5],
         indiv.tab[seq(2,15,3), 2],
         indiv.tab[seq(2,15,3), 5],
         indiv.tab[seq(3,15,3), 5],
         indiv.tab[seq(3,15,3), 2])

up <- c(indiv.tab[seq(1,15,3), 3],
        indiv.tab[seq(1,15,3), 6],
        indiv.tab[seq(2,15,3), 3],
        indiv.tab[seq(2,15,3), 6],
        indiv.tab[seq(3,15,3), 6],
        indiv.tab[seq(3,15,3), 3])

segments(x.coords, low, x.coords, up)
points(x.coords, est, pch=(rep(c(2, 1, 16, 0, 15), 6)), col="black", bg="black", cex=1)
abline(h=0.653, lty="dashed")

mtext("Gender", outer=TRUE, cex=1.3, at=0.20, line=-1, font=2)
mtext("Party", outer=TRUE, cex=1.3, at=0.5, line=-1, font=2)
mtext("Education", outer=TRUE, cex=1.3, at=0.835, line=-1, font=2)

mtext("Male", outer=TRUE, cex=1, at=0.13, line=a, font=2)
mtext("Female", outer=TRUE, cex=1, at=0.27, line=a, font=2)

mtext("Republican", outer=TRUE, cex=1, at=0.44, line=a, font=2)
mtext("Democrat", outer=TRUE, cex=1, at=0.58, line=a, font=2)

mtext("No Higher", outer=TRUE, cex=1, at=0.76, line=a, font=2)
mtext("Higher", outer=TRUE, cex=1, at=0.9, line=a, font=2)

legend("bottomright", pch = c(2, 1, 16, 0, 15), 
       legend = c("Direct Questioning", "List Experiment", "List Experiment with Auxiliary Information", 
                  "Randomized Response", "Randomized Response with Auxiliary Information"), 
       cex = x, 
       bty = "n")
dev.off()

